Tests of Dynamical Scaling in 3-D Spinodal Decomposition 
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We simulate late-stage coarsening of a 3-D symmetric binary fluid. With reduced units l,t (with 
scales set by viscosity, density and surface tension) our data extends two decades in t beyond earlier 
work. Across at least four decades, our own and others' individual datasets (< f decade each) 
show viscous hydrodynamic scaling (I ~ a + bi), but b is not constant between runs as this scaling 
demands. This betrays either the unexpected intrusion of a discretization (or molecular) lengthscale, 
or an exceptionally slow crossover between viscous and inertial regimes. 
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When an incompressible binary fluid mixture is 
quenched far below its spinodal temperature, it will phase 
separate into domains of different composition. For sym- 
metric (or nearly symmetric) mixtures, these domains 
will, at late times, form a bicontinuous structure, with 
sharp, well-developed interfaces. The late-time evolu- 
tion of this structure in three dimensions remains incom- 
pletely understood despite theoretical (l]-|j| experimental 
H and simulation |5]H work over recent years. 

In the present work, we use the DPD (dissipative par- 
ticle dynamics) simulation algorithm || to access length 
and time-scales far beyond those reached previously. (De- 
tails of the simulations will appear elsewhere |fl0|| .) When 
combined with other datasets |-|^] our results allow a 
severe test of the dynamical scaling ideas which underlie 
most theoretical treatments jl]-3| and data analyses |Q. 
We conclude that dynamical scaling is in doubt, perhaps 
due to the intrusion of a molecular lengthscale through 
the physics of topological reconnection events. An alter- 
native explanation of the results, based on a universal 
but extremely slow crossover, is also carefully examined. 

As emphasized by Siggia jij, the physics of spinodal 
decomposition involves capillary forces, viscous dissipa- 
tion, and fluid inertia. Indeed, assuming that no other 
physics enters, then the parameters governing the behav- 
ior are the interfacial tension a, fluid mass density p, and 
viscosity rj. (We now specialize to 50/50 mixtures with 
complete symmetry of the two species. Any asymmetries 
in composition, thermodynamics or viscosity [[tl| provide 
additional control parameters.) From these three param- 
eters can be constructed only one length, L = r] 2 /pa and 
one time T — rj 3 /pa 2 . We now define the lengthscale 
L(T) of the domain structure at time T via the structure 
factor S(k) as @ L = (2tt) (/ kS(k)dk/ j S*(fc)dfc) _1 . 
The exclusion of other physics in late-stage growth then 
leads us to the dynamical scaling hypothesis 
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where we define reduced time and length variables via I = 
L/Lo and t = T/Tq. Since dynamical scaling should hold 
only after interfaces have become sharp, and transport 
by molecular diffusion suppressed, we have allowed for 



a nonuniversal offset a in eq.[l]. Thereafter the scaling 
function f(t) should approach a universal form, the same 
for all (fully symmetric, deep-quenched, incompressible) 
binary fluid mixtures. 

It was argued further by Siggia jlj that, for small 
enough t, fluid inertia is negligible compared to viscos- 
ity, whereas for large enough t the reverse is true. This 
imparts the following asymptotes to the function /: 
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where, if dynamical scaling holds, the amplitudes b and c 
must be universal, as must the crossover time t* (defined 
e.g. by the intersection of asymptotes on a log-log plot). 
Note that the Reynolds number Re = (pL/r))dL/dT = 
ff which becomes large in the inertial regime, eq.|^. 

Perfectly symmetrical fluid pairs do not exist in the 
laboratory, but computer simulations allow us to test the 
validity of eq.l, on which the wider interpretation of ex- 
periments crucially depends In the viscous regime 
(eqj|), the scaling reduces to L(T) = A + BT where A is 
nonuniversal and B = ba/ri- This linear law has been re- 



ported by several groups J 
only in two recent cases t 



0,| (see also HIH) but 
were reliable a and rj values 
obtained, as are needed to find b. In both of these, the 
offset A was significant, and the linear regime (straight 
part of the curve at late times) spanned much less than 
a decade. In reduced units |l^] , we find that the data of 
Ref. describes times in the range 1 < t < 3 with a 
value of b = 0.3. However the MD data of Laradji et al 
§ has 60 < t < 140 and b = 0.13. 

The discrepancy over b (see also [^)) cannot simply be 
brushed aside. For if dynamical scaling (eq.^j) applies, 
and both simulations [gjj^] are (as claimed) in the viscous 
regime (eq^), then these two b values should both be 
the same |[L7 |. It is thus premature to conclude that any 
universal regime of viscous hydrodynamic scaling (eq.|^) 
had yet been observed in computer simulations. 

To clarify this important issue, we have conducted sev- 
eral new simulations which vastly extend the range of 
timescales explored: we probe 750 < t < 45, 000. This 
was done using the DPD algorithm which combines soft 
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interparticle repulsions with pairwise damping of inter- 
particle velocities and pairwise random forces ||. The 
latter conserve momentum, leading to a faithful simula- 
tion of the isothermal (and, in this study, effectively in- 
compressible JlTl ) Navier Stokes equation at large length- 
scales. Among several advantages of DPD over MD, ex- 
ploited below, is that the viscosity of a DPD fluid can be 
varied independently of its thermodynamics. 

To describe our DPD parameters, we briefly switch 
from reduced physical units (l,t) to "DPD units": the 
range of the repulsive interaction is unity, as is the par- 
ticle mass. We further set ksT = 1 (T is temperature). 
With the form of repulsion used by Groot and Warren 
||, we chose a particle density 10 and energy parame- 
ters an = &22 = 20, ai2 = 100, which is a deep quench 
(T c /T ~ 80) Jl(J. The timestep was 0.01 §, giving mea- 
sured T's within 2 percent of the nominal value. Inte- 
grating the microscopic stress across a flat fluid-fluid in- 
terface the interfacial tension was found as 50.6±0.2. 
For each damping parameter 7 the viscosity was found 
from the the mean stress in steady shear under Lees- 
Edwards boundary conditions Jl9[ ] ; values varied between 
77 = 2.6 ± 0.2 (7 = 1) and 77 = 12.2 ± 0.5 (7 = 30) g(J. 

Most runs were performed on a 512 node Cray T3D, 
with a typical run time of several thousand processor- 
hours. Resources allowed one or two full-sized runs for 
each viscosity. The simulation box for these contained 
10 6 particles with periodic boundary conditions. Thor- 
ough tests of scaling and data collapse for S(k) were 
made (l(J. Finite size effects became apparent when the 
structural lengthscale L exceeded about half the box size 
(L > A/2 ~ 20); data beyond this was excluded from 
our fits for f(t). We also excluded an "early stage" por- 
tion of each run; this was judged by eye from the shape 
of the L{T) plot. (Possible resulting bias is considered 
below; little would be changed had we instead applied a 
sharpness criterion to the observed interfaces.) 

The datasets for L(T) (DPD units) are presented in 
Fig.l (inset). Excluded early time data is shown dot- 
ted, as is some data for L > A/2. Slight wobbles in the 
fitted parts of the curves represent sampling errors in 
L arising because L/A is not small; these vary between 
duplicate runs and appear distinct from the direct finite 
size (saturation) effects arising for L > A/2 [^,|T). Fig.l 
shows the same data (with offsets subtracted) on a log- 
log plot. Note first that, since most of the plots in Fig.l 
show upward curvature at early times, our elimination of 
early time data will bias downward any estimate of the 
quantity z — d In f /d In t (a true or effective scaling ex- 
ponent). Despite this, only for the smallest viscosity run 
(if that) is there appreciable direct evidence for an expo- 
nent z < 1 as predicted for t ^> t* (eq.||). A fit to eq.|l| 
with / = ct z in fact gives z = 0.88 whereas all but one 
of the other viscosities give 1.10 < z < 1.17 (77 = 9.8 has 
z = 0.96). This suggests that our lowest viscosity run 
(r/ = 2.6) and it alone, may be approaching the inertial 
hydrodynamic regime, eq||; for more evidence of this see 
0. This run covers 20, 000 < t < 45, 000, implying that 
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FIG. 1. Inset: Raw DPD data; L vs. T for viscosities (left 
to right) n = 2.6,3.5,4.6,6.2,8.2,9.8,12.2. The datasets for 
77 = 6.2, 9.8 are averages of two runs. Main figure: the same 
data in reduced units (log-log) with offsets (found by linear 
extrapolation to t = 0) removed. 



t* (the crossover between eqsj2|j^) is similarly large. A 
less extreme number is obtained if one quotes instead the 
equivalent Reynolds number Re* ~ b 2 t* . (This relation 
applies because in linear scaling, we have Re = bl ~ b 2 t.) 
For the middle of the given run, Re is about 20, so Re* 
need not be much larger than this. Given the smallness of 
the apparent b values (see below), the largeness of t* fol- 
lows, as does the failure to observe a clear inertial scaling 
regime (eq.||) in previous simulations ||. 

Based on these observations, we have fit our remain- 
ing 6 datasets to the viscous hydrodynamic scaling form, 
eq.|[ In all cases the fits are at least as convincing as 
those of 0,|t| . Despite this, we definitely cannot interpret 
this data (nor that of |^,0|) as support for a universal vis- 
cous hydrodynamic scaling, eq.|^. Fig. 2 (a) shows fits to 
/ = bt (deviations from the data invisible on this scale) 
including our own and earlier ||[7j datasets. Fig. 2(b) 
shows the fitted b coefficients against the mean time i, 
defined by the middle of the fitted section of each run. 
Obviously, b is not constant as required: it drifts system- 
atically toward smaller values at later times t, a trend 
represcntablc empirically as a weak power law, b ~ f~ 02 . 

What does all this mean? Clearly, one would expect 
to measure b ~ i~°- 2 if in fact one had / = ct z with 
z ~ 0.8 (see Fig. 2(a)). We are aware of no theory pre- 
dicting this value of z, so it would presumably have to be 
interpreted as an intermediate, effective exponent arising 
in the crossover region between eqsj^ and [|. Although 
possible, at least two arguments counter this interpreta- 
tion. Firstly, the "crossover" , if this is indeed what we 
are seeing, must be exceptionally broad. Fig. 2(b) shows 
that a single effective exponent governs the entire range 
of data shown: any "crossover" region covers four decades 
in time. The second reason to doubt this explanation is 
that for all of the DPD datasets shown in Fig. 2, a fit to 
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/ = ct z yields values of z that are not close to 0.8, but 
close to (and usually slightly larger than) 1.0. Put differ- 
ently, even after subtraction of the offsets a, our datasets 
do not join up into a continous curve on the (I, t) plot as 
dynamical scaling requires. This is apparent from Fig.l, 
and remains true under various replottings we have tried 
(such as recalculating offsets by imposing z — 0.8 rather 
than z = 1). We therefore ask whether there might be 
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FIG. 2. (a) Fitted functions f(t) = bt for DPD data 
(rightmost 6 datasets) that of Ref. [6] (centre left) and Ref. 
[7] (far left), (b) log- log plot of resulting growth velocities b 
against the midpoint time t of each run. 



some other physics, playing a role in spinodal decompo- 
sition at late times, which could lead to a violation of 
the dynamical scaling hypothesis itself. One possibility 
is that the late-stage coarsening velocity b depends on 
initial conditions, inherited from the nonuniversal early- 
stage dynamics. (For related ideas, see plfl.) This infor- 
mation would have to reside either in the velocity field 
itself, or in subtle details of the density distribution. The 
first of these can be tested numerically by re-initializing 
the fluid velocity field during a late-stage run; we have 
done this and no significant effect on b was observed. 

A more plausible mechanism for the observed nonuni- 
versality of the velocity b could arise from the direct in- 
trusion of physics that the dynamical scaling hypothesis 
excludes. Thermodynamics (e.g. finite temperature or 
compressibility fL7| ) cannot be solely responsible, since 
all our DPD runs are identical thermodynamically. Per- 
haps the most interesting possibility is that late-stage 



spinodal decomposition involves a molecular (or, in simu- 
lations, discretization) lengthscale which could enter dur- 
ing topological reconnection or "pinch-off" events. In 
such events, without which coarsening of a bicontinu- 
ous structure cannot proceed, a fluid neck contracts to 
(formally) zero width in finite time. 

Recent work on a closely related problem (disconnec- 
tion of a single fluid domain in vacuo) suggests that 
pinch-off processes need not violate dynamical scaling 
fl23fl : the asymptotic behavior both before and after 
the pinch have a universal description in l,t variables 
(measured from the pinch-off event itself). According 
to this work, molecular physics intervenes only briefly 
at pinch-off, and is forgotten soon after. It is not yet 
known whether similar universality can be recovered for 
fluid- fluid pinch-offs [^3| , but crucially, even in the fluid- 
vacuum case, such universality is only expected for large 
values of the dimcnsionless quantity 



A = Lq/H — rj 2 / pah 



(4) 



where h is a molecular (or discretization) length [p3|,p[. 
For the fluid/vacuum case, Eggers J2i| argues that A is 
large enough for some fluids (~ 10 7 for glycerol) but not 
others (~ 20 for water), to recover universal behavior. 

If similar ideas govern the fluid-fluid case, and if pinch- 
off physics remains a controlling factor in late stage coars- 
ening, then a violation of dynamical scaling could be ex- 
pected for many real fluids. The same applies for any sim- 
ulation in which A is not very large. Taking h = 1 (DPD 
units) we find that A in our runs ranges from A = 0.28 
at r\ = 12.2 (so that t = 800) to A = 0.014 at r) = 2.6 
(so that i = 30000). The systematic dependence of b on 
t reported above can, for these DPD runs, equally well 
be expressed as a dependence on A. The latter would 
permit an extended form of dynamical scaling, with f(t) 
replaced by /(A, t) in eq.0; at present this cannot be dis- 
tinguished from a t dependence, because the variations 
we make through rj affect i and A similarly [p4| . 

One speculative possibility is that the time AT taken 
for a fluid neck, of order the domain size L, to reach 
pinch-off is not linear in L (as eq.|| suggests) but varies 
as L\n(L/h) (so, in reduced units, At ~ Zln(ZA)). In this 
case, individual runs would show little departure from 
f = bt, yet b would drift slowly downward with i, and 
runs of different rj would not quite superpose on the (I, t) 
plot. Such logarithms could conceivably arise from the 
hydrodynamics of thin fluid cylinders |25(| . A fit to b — 
ln?7 (not shown) is comparable in quality, for our DPD 
runs, to that in Fig. 2(b) but less good than the power 
law shown there, if the data of |||7]] is included. 

In conclusion, we have made a careful analysis of our 
extensive new DPD data Jl0[ | and of previous simulation 
results on spinodal decomposition in fully symmetric 
binary fluids. Taken together, these data now cover ap- 
proximately five decades in reduced physical time units. 
Contrary to expectation, the data offer no clear support 
for the hypothesis of a universal dynamical scaling (eq.[l]) 
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pq . Such a hypothesis can be sustained, but only M] by 
assuming an extremely broad crossover between viscous 
(z = 1) and inertial (z = 2/3) scaling regimes, with an 
effective exponent z ~ 0.8 spanning about four decades 
of reduced time t. This "slow crossover" interpretation 
is more plausible when expressed in terms of Reynolds 
number, which spans only 0.1 <Re< 20; indeed, expe- 
rience with turbulence |[27f shows that a wide regime of 
Re might exist in which inertial effects are significant 
(spoiling eq.||) but not dominant (as required for eq.^J). 
However, such an interpretation leaves unexplained the 
facts that (i) almost all our individual simulation runs 
are better fit by z > 1 than z ~ 0.8, and (ii) even after 
subtraction of nonuniversal offsets, these runs do not lie 
on a common curve on the I , t plot |2l|] . 

We have therefore put forward a more radical proposal: 
that the list of physical ingredients (fluid viscosity, den- 
sity and interfacial tension), assumed by the dynamical 
scaling hypothesis to dominate the physics of spinodal 
decomposition at late times, is incomplete. One possibil- 
ity is that the ratio A = Lo/h (of continuum to micro- 
scopic lengths) remains a relevant parameter: coarsening 
of a bicontinuous structure is contingent on topological 
reconnection (pinchoff) events, which could allow the in- 
trusion of microscopic physics no matter how large the 
mean domain size L. 
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